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Abstract 

In the field of Large Eddy Simulation, the Smagorinsky subgrid scale model (in some 
form) is the most commonly accepted and used subgrid scale model. The purpose of 
this paper is to address the main weakness of the Smagorinsky model, its poor perfor- 
mance near the wall. The goal is to establish a model that corrects the Smagorinsky 
model near the walls while at the same time minimizing the computational over- 
head. A version the Dynamic Subgrid Scale model is also incorporated into the finite 
element code to facilitate comparisons with the new model near the walls. 

One of the unique characteristics of Large Eddy Simulations as compared to other 
methods of dealing with turbulent flows is the idea of filtering. In this paper we define 
what a filter is and also address an issue related to filters; the error that results when 
the filtering and differential operations are interchanged. This error is studied under 
the context of the Finite Element Method which allows us to focus on the function 
being filtered rather than the filter kernal function, which has been the usual approach 
in studying this error. 
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1 Introduction 

In the study of turbulence, three distinct computational approaches have been devel- 
oped to simulate turbulent flows. Direct numerical simulation (DNS), as the name 
implies, is the most straight forward approach to the simulation of turbulent flows. 
The DNS approach is to simply solve the Navier-Stokes equations using some nu- 
merical method. The problem with this approach is that to succesfuUy solve the 
Navier-Stokes equations, the computational domain must be large enough to contain 
the largest scales of motion, L, while the grid resolution must also be fine enough 
to resolve the smallest scales, which are of the order of the Kolmogorov microscale, 
?7 = (i/^/e)-*^/^, where rj is the dissipation rate per unit mass and v is the kinematic 
viscosity. So, the number of grid points necessary in each direction can be estimated 
as ^ and using the well established relationship between ^ and the Reynolds number. 
Re, (see Tennekes and Lumely, [26]) we have that 



where N is the number of grid points. Hence, it is clear that DNS is only practical 
for low Reynolds number flows. At the other end of the spectrum from the DNS 
approach is the Reynolds averaged simulations (RAS). In RAS, the flow (velocity fleld 
and pressure for incompressible flows) is broken up into a statistically steady (mean) 
portion and fluctuations. The mean flow is then solved for, while the effects of the 
fluctuations on the flow are modelled. By modeling the fluctuation effects, the RAS 
approach models all the turbulent interactions whereas the DNS approach completely 
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resolves them. One of the main stumbhng blocks for RAS is that no one model has 
been found for the fluctuation effects that can be used for all turbulent flows. It seems 
unlikely that such a model would exist given the complexity of turbulent flows and the 
fact that turbulence is a property of the flow rather than the fluid. The third approach 
to the computation of turbulent flows, known as Large-Eddy simulation (LES), is 
often viewed as the intermediary approach between DNS and RAS. The motivation 
behind LES is that since the large energy carrying eddies are highly influenced by the 
boundary conditions, it should be computationally resolved, while the small eddies 
or unresolved scales are modelled. It is hoped that since the small eddies are more 
homogeneous and isotropic, a simpler and more universal model can be used for the 
unresolved scales as compared to the RAS models. 

The objective of this study is to develope a new model for the unresolved scales 
that overcomes some of the problems associated with the established models, which 
are discussed in the coming sections. This study is based on and, in many ways, 
an extension of the work begun by Rose McCallen [19] in which a LES method 
is incorporated into a finite element code. In her study, a classical model for the 
unresolved scales developed by Smagorinsky [25] was used to simulate the channel 
and backward facing step flows in two dimensions. In the present study, this model 
plus a version of a newer model first developed by Germano, Piomelli, Moin, and 
Cabbot [6] and our own model base on the Smagorinsky model are used to simulate 
the channel and backward facing step flows in three dimensions. 
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2 LES 

2.1 Filtering 

In LES, the underlying principle is to seperate the large eddies, which are computa- 
tionally resolved, from the small, which are modelled. This partitioning of the scales 
is accomplished by filtering, in the case of incompressible flow, the velocity field, u, 
and pressure, P. Following the work of Leonard [16], if /(x) is a function that con- 
tains all the scales, then the filter of /, denoted as /, is defined as the convolution of 
/ with a filter function G{x). 

7(x) = I^G(^-s)f(s)ds (1) 
where F is the flow volume, G is normalized so that 

J^G{s)ds^l, 

and G{—x) = G{x). This is necessary to insure that when / is constant that / = /. 
Some common filters are the following: 



• Top-Hat filter 

ior|x — s| < ^ 

^(x - s) = < 



^ for|x-s|<^ 



otherwise 

where A/ is the filter size. If u is the velocity field, then the filter of u, using 
the Top- Hat filter, is 

i / 2 



u(x) = ^ y^I^ u(x + s) 
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which is equivalent to cell volume averaging of Deardorff [2] . Now, the Fourier 
transform of u(x) is 



where u(k) is the Fourier transform of u (see Kwak et al. [14].) So, we see 
that the spectrum of the filtered field contains components of all wave numbers. 
This implies that the Top-Hat filter will not filter out all the small scales of 
the flow. Also, the wave numbers for which the coefficient of u(k) in Eqn. (2) 
is zero lead to the inverse transform being singular. This means that it is not 
possible to obtain the actual spectrum u(k) from the filtered one, u(k). So the 
Top-Hat filter is not appropriate when spectral results are sought. 

• Gaussian filter 



According to Rogallo and Moin [24], the Gaussian filter is a good filter choice for 
filtering in the homogeneous directions because it provides a smooth transition 
between the resolved and subgrid (unresolved) scales and is positive definite in 
both the physcial and wave space. Kwak et al. [14] showed that the Gaussian 
filter results in a filtered field that captures most of the large scale motion while 
most of the small scale motions are removed. This ability of the Gaussain filter 
to seperate the large and small scales, gives it a desirable property for a filter 
in LES. The one obvious fiaw or weakness of this filter is that it does not have 
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a compact domain. Hence, its use in non-homogeneous, wall bounded flows is 
questionable at best. 

• Sharp-Cutoff filter 



Just as the Top-Hat filter has compact support in grid coordinates, the Sharp- 
Cutoff filter has compact support in the wave space. Also, just as the Top-Hat 
filter looses it compact support when transformed into the wave space, the 
Sharp-Cutoff filter looses its compact support when going from the wave space 
to grid coordinates. So, like the Gaussian filter, it is not appropriate for non- 
homogeneous, wall bounded flows. 

Although Leonard's definition for filtering is generally accepted, two distinct ap- 
proaches to the apphcation of the filter have been developed. These are mainly 
infiuenced by the fiow under study and the numerical technique used to solve the 
flow problem. In the case of homogenous turbulence, the use of spectral methods 
facilitates the use of an explicit filter. However, in using finite difference methods, 
usually for non-homogeneous fiows, the filter is assumed to be implicitly defined by 
the discretization. An example of this implicit filter is given by Rogallo and Moin, 
[24] , where they note that the second-order central-difference of a function / is the 
derivative of the filter of /: 
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where the filter used is the Top-Hat filter. Some of the reasons for using an implicit 
filter rather than explicitly applying a filter are the following: the difficulty of actually 
applying a filter; the reduction in grid resolution resulting from the fact that in 
finite differencing, the function is only defined at the grid nodal points and so any 
filtering would require at least two nodal points (i. e. just like the pressure, the 
differencing would be done using the grid nodes, but the filtered value would be at the 
cell center) ; and finally, finite difference schemes are usually used for non-homogenous 
wall bounded fiows which usually have anisotrophic grids (since the filter width will 
most likely depend on the grid, this means that the filter width would be a function of 
location unlike the homogeneous case where the filter width remains constant. This 
leads to further difficulties, as will be discussed below.) 

In applying a filter, Eqn. (1), to get the LES equations, a question that arises is 
whether or not the derivative operator and filtering are interchangeable, i. e. is the 
derivative of the filter equal to the filter of the derivative. The answer will depend 
on the flow under study. In the case of homogeneous turbulence, where the scales of 
the large and small eddies do not depend on the location in space, there would be no 
reason to use anything but a constant filter width, A. So, assuming that G ^ as 
t ^ oo and using the 1-dimensional case for notational simplicity (with the extension 
into higher dimensions following in the usual way). 
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where integration by parts and the fact that ii ^ — x — s, then 



dx dx 



and 




were used in the above derivation. So, under the assumptions made above, the deriva- 
tive and filtering operations are interchangeable in the case of uniform filter width 
which occurs in homogeneous turbulence. 

For the non-homogeneous wall bounded flow, since the scale of the large and small 
eddies will depend on position relative to the wall, it is clear that the filter width must 
be a function of location. To emphasize this, we rewrite the definition of filtering, 
Eqn. (1), as 



It can be shown by simply carrying out the differentiation, that the derivative and 
filtering operations are not, in general, interchangeable. For this reason, there is 
ongoing research to find filters that reduce the error in interchanging the derivative 
and filtering operations (Germano [5]; Ghosal and Moin [7]; van der Ven [28]; Najjar 
and Tafti [21]; Vasilyev and Lund [30] [29]; and Marsden, Vasilyev and Moin [17].) 
However, as of the time of this report, there is no compelling evidence to indicate 




(3) 
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the best filter choice. As mentioned before, in using the finite difference methods, 



researchers often assumed an implicit filter. The issue of the interchange of filter and 
differentiation was then bypassed by assuming that just as when using the Reynolds 
Average method, one assumes that all turbulent behavior is captured by the model 
for the Reynolds stress, all of the effects of filtering were assumed to be captured 
by the sub-grid scale model. More recently, there have been some studies where an 
explicit discrete filter is used instead of the implicit filter (Najjar et al. [21]; Vasilyev 
et al. [30] [29]; and Marsden et al. [17].) However, as Vasilyev pointed out, there is 
more work to be done before any conclusions can be drawn about the effectiveness of 



the discrete filter. 



In the case of the Top-Hat filter, in 1-dimension, the filter of / is 

1 



/ 



AM 
2 



A{x) Jx-^ 



f{s) ds 



and 



1 



_ I A(£) 
2 



df 



dx A(x) Jx-^ 



] ns)ds 



A(x) 



Using the above, we calculate the derivative of the filter of / as follows: 



dl ^ d_ 
dx dx 



A{x) Jx 

A(£) 

\' f{s)ds + 



2 

A'(x) . . 1 



A2(x) .X 

1 



Me) 

2 



A{x) 



A(a;; 
A(x 



df A'jx) A'jx) 
dx A{x)-' 2A{x) 
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By the trapezoidal rule, 



^ A(£) ^ 



f{x + 



A(x) 



+ f{x 



A(x) 



+ 0(AY') 



So, using this and the previous equation, we get 



Since we choose our filter width. A, to be the length of our grid element and the 
grids that we use are graded near the wall. A' ^ and so the above shows that the 
magnatude of the error in the interchange of the derivative and filtering operator is 
dependent on the derivative of the filter width. A'. To obtain a better prespective on 
the implications of the interchage error, we first rewrite it as 



dx dx ^A ■' ' 



It is clear that if ^ = then the interchage error is second order. Now, assuming 
that 



A^ 

'a 



= M 



for some constant, M, we get 



A = Ce 



Mx 



where C is a constant. Hence, if the filter width behaves like the exponential function, 
then the error in the interchange of filtering and differentiation is second order, using 
the Top-Hat filter. 
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2.2 Filtering and FEM 

One of the advantages and/or characteristics of the Finite Element Method (FEM) 
over the Finite Difference Method is that instead of producing a solution over a 
discrete set of points, it produces a continuous approximation of the actual solution. 
In fact the general shape of the approximation function is known and chosen by the 
user via the shape functions before the actual computations begin. This knowledge of 
the general shape of the function being filtered naturally leads to the question of how 
this information can be used or effects the analysis of the error in the interchange of 
filtering and differentiation. 

We begin by modifying the definition of filtering, for the one dimensional case, in 
the following way. 

1 f^+^ fx-s\ 
u(x) — ^ , , / , G\ Ue(s) ds 

where Ue{x) is the finite element representation of u{x). Now using the transformation 
y — the above becomes 

1 

u{x) ^ \ G{y)ue{x - yA{x)) dy 

Differentiating this, we get 

^(^) = G{y)u'Xx - y/\{x)){l - yA'{x)) dy. (4) 

Note that in general Ue{x) is only required to be continuous rather than continuously 
diffcrcntiablc. It is however continuously differentiable in each element. So we use 
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the above notation with the understanding that 



f\ G{yK{x - yA{x)){l - yA'{x)) dy = T G{yK^{x - yA{x)){l - yA\x)) dy 

/ G{y)u',^{x-yA{x)){l-yA'{x))dy 



+ 



+ 



+ 1'' G{y)v!,^{x-yA{x)){\-yN{x))dy 



where (x^, Xi+i) represents the transformed elements contained in the support for the 
filter. Note that 

^{x) = J^^G{yH{x-yA{x))dy 
and so using this in Eqn. (4) we get our error term 



du du 



A\x) / ^yG{y)u'^{x - yA{x)) dy. 



(5) 



If we assume a linear shape function then 

Uei{x) = ai + Pi{x) 

for X E where is the i-th element. We define the filter width A{x) to be twice 
the width of the smaller of the two intervals containing the nodal point x. This allows 
us to seperate the error term, Eqn. (5), in terms of each element as follows: 



du , , du , , 



(3,A'{x) yG{y)dy + (52A\x) J_^yG{y)dy. 



The above would seem to indicate that the error in the interchange is of order A'{x). 
This would be the case if the filter width is defined such that a change in the filter 
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width does not imply a change in the element width. However, the above definition 
of the filter width not only involves x but also the element width, which we denote 
Af,i{x). Now, we note that 

ft = = «'w +o(A„(x)) 



and 



P^ = u'{x)+0{A,,{x)) 



where if we let x = Xn, then (x) = Xn — Xn-i and A^^ (x) = Xn+i — x^ for gradings 
of Xn that increase with n. Since G is defined as an even function we have 

yG{y) dy = - I yG{y) dy. 



So 



^(^)-^(^) = W2-^)A\x)j^'yG{y)dy 

(0(Ae,(x)) - 0{A,,{x)))A'{x) r yG{y) dy 

Jo 

= 0{Aeix)A'{x)) 

where for practical problems the order of the local element widths are the same, 
i. e. 0(Aei(a;)) = 0{Ag2{x)) = 0{Af,{x)). Furthermore, if we assume a smooth 
exponential grading such that 0{A'{x)) — 0{A{x)), then the interchange error is 
second order, 

|M-fw^0(A^(.)). 
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For a quadratic shape function, we have 

Uei{x) = Ctj + PiX + Jix"^ 

for X e Again, using the above filter definition, the interchange error is 



du du 



(/3i + 27ix)A'(x) ryG{y)dy 
Jo 

27iA(x)A'(x) ry'G{y)dy 
Jo 

{(32 + 2j2x)A'{x) [\yG{y)dy 
272A(a;)A'(a;) [\y'G{y)dy 



Note that in FEM, Uei{x) = ai + PiX + jix'^ is exact at each nodal point, i. e. Uei{xn) ~ 
u{xn)- So, using the above definition of Ae;(a:), we have the following first order 
approximation of the derivative of u: 

, u(x) — u(x — Ae,(x)) ^, . , 

^'(^) = ^ ^ ^^^^^ ^^ +0(Ae,(x)) 

Ue^jx) - u^^jx - A^^{x)) ^ 
_ Oil + Pix + ^ix^ - [ai + Pi{x - Aej^jx)) + 7i(a: - Aei(a:))^ ^ 

Aei(x) 

= A + 271a; - 71 Aei (X) + 0( Aei (X)) 

So 

/3i + 27ix = M'(x) + 0(Aei(x)) 
and using a similiar argument, we also have 

/32 + 2-f2X = u'{x) + 0{A,,{x)). 



+ 0(Ae,(x)) 
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Hence, using the symmetry of G as in the hnear case, the interchange error is again 



du . . du . . 



0(Ae(x)A'(x)). 



We now show that for any polynomial shape function, the interchange error is of 
order 0(Ae(x) A'(x)). Consider a general n-th order polynomial shape function 



i=0 



Note that 



u 



i-l 



i=l 



So the interchange error, Eqn. (5), can be written as 



du du 



+ 



A'{x) lyG{y)J2tai,,{x-yA{x)y-Uy 

•'~2 i=l 

A\x) r yG{y) ^«.2(x - yA{x)y-' dy 



We can ignore any term that has A{x) since they will be of oder 0{A{x)A'{x)) 
or higher. So, using the Binomial Theorem and the fact that G is symmetric, the 
interchange error above can be rewritten as 



du . , du . , 

Tx^""^ - Tx^""^ 



A\x)[{Y^ia,,,x'-^) - {Y.^a,,i^"')\ lyG{y)dy (6) 



i=l 



i=l 



+ 0{A{x)A'{x)) 



Now, as in the quadratic case, we note that the interpolating polynomial function is 
exact at each nodal point, i. e. Ue.{xn) — u{xn). So 

u{x) — u{x — Aei(a;)) 



u'{x)+0{Ae,{x)) = 



Aei(x) 
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Ae,{x) 

Aei(x) 

Again, using the Binomial Theorem, we rewrite the above as 

, , E7=o«M^^-Er=o«.,iE}=oC>^(-Ae,(x)r^' 
u [x) + 0{Ae,{x)) = ^ 

Note that the terms that do not have any factor of Ag^{x) cancel out which leaves 
the terms with a factor of just Aei(a;) as the lowest order terms. So, collecting these 
terms, the above can be rewritten as 

n 

u'{x) + OiA,,{x)) =J2m,ix'-\ 

i=l 

Using same argument as above only with a first order forward differencing scheme 
instead of the backward differencing scheme, we get 

n 

u'{x) + 0{A,,{x)) =J2m,2x'-'. 

i=l 

Using the above in Eqn. (6), we get 

= 0{Ae{x)A'{x)) 



du du 



where again we assume that the filter and element widths are of the same order, i. e. 
0{A{x)) = 0{Ae^{x)) = 0{Ae{x)). So if we assume that the grading of the elements 
is smooth enough such that 0{A'{x)) = 0{A{x)) then the error in the interchange of 
differentiation and filtering is second order for any polynomial interpolation of degree 
one or greater. 
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Appendix A presents the results of a computational experiment that verifies the 
above results for the case of linear and quadratic shape functions. Note that since 
the FEM is apphed to solve the filtered equations, the above filtering technically only 
applies to filtering a filtered quantity. This is true for any explicit filter of u since 
the only values of u that are available are those that are obtained from the numercial 
method. But the numerical method is used to solve the LES or filtered equations 
and not the original, as discussed in the following section. Hence, to be completely 
consistent, an explicit filter can only be applied in a double filtering situation such 
as when using the Dynamic Subgrid Scale model discussed below. It would seem 
that to be consistent the initial filtering needs to be introduced in one of two ways. 
One, it can be introduced by a term representing the derivative/filter interchange 
error or this error can be ignored by using a filter whose interchange error is of the 
same order as the numerical scheme. Note that this error must be derived using the 
original vector field u and should avoid any reference to the computed u. Also, this 
filter would not be applied directly to the numerical solution but should instead be 
applied to the experimental or DNS data for comparison with the LES data. The 
second and ideal place to introduce the filter effect would be the SGS model since the 
filter determines the scales that are in the system. 
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2.3 Governing Equations 

The set of equations that govern the motion of (Newtonian) fluids are known col- 
lectively as the Navier-Stokes equations. In the case of incompressible fluids with 
uniform density, the non-dimensionalized equations are 



where we use the standard tensor notation of repeated indices indicating summation, 
P — J, and Re—^ with U begin the characteristic velocity, L begin the character- 
istic length, and u being the kinematic viscosity. Now, assuming that the error in the 
interchange of filtering and differentiation can be ignored (i. c. |^ = |^), the filtered 
Navier-Stokes equations are 



By filtering the equations above, we have introduced a new unknown term, UaUp, but 
no new equation. This means that the system of equations is now underdetermined. 
This is known as the closure problem, which is also encountered by researchers using 
RAS. In fact, the Reynolds Averaged Navier-Stokes equations look very much hke 
the filtered Navier-Stokes equation above (Eqn. (7) with the " replaced by < • >, 
which represents the Reynolds averaging operator.) The difference is in how the " 
(and/or < • >) is defined. In LES, the decomposition of variables is in terms of the 
large scales and small, i. e. u = u + u' where u is the large scale component of u and 
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u' is the small or sub-scale component. Prom the definition of filtering, Eqn. (1), 
it is clear that in general u ^ u, unless of course u is constant or the filter used is 
the Sharp-Cutoff. This means that, in general, u' 0. In the Reynolds Averaging 
method, the decomposition is in terms of the statistical mean and fiuctuations with 
mean. So if we write u =< u > +u' then << u »=< u > and < u' >= 0. This 
difference, in the interpretation of the decomposition of variables, leads to a very 
different interpretation, both physically and mathematically, of Eqn. (7). In either 
case, note that the interaction between large and small scales {u and u') or mean and 
fiuctuations {< u > and u') must occur in the u^/up term of Eqn. (7). As this work is 
based on LES, our focus will be on the closure problem as it relates to LES. 

To deal with the closure problem, the above filtered equations, Eqn. (7), are 
rewritten as 

1^ = 

(8) 



dt ^ dxg^ °' dxc ^ ReS^^"° 



where 



is called the subgrid scale (SGS) Reynolds stress. The SGS Reynolds stress in LES is 
simihar to the Reynolds stress in RAS in that in both cases the the unresolved scales, 
u' (small eddies in LES and turbulent fiuctuations in RAS,) are viewed as producing 
stresses in the resolved scales (large eddies in LES and the mean velocity in RAS.) 
The difference is that in RAS, u' represents all the turbulent motions, while in LES, u' 
represents only the small eddies or subgrid scales. This means that the energy in the 
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unresolved scales of LES (subgrid scales) is a much smaller portion of the total flow 
compared to the energy in the unresolved scales of RAS (the turbulent fluctuations.) 
Hence, it is believed that the accuracy in the modelling of the SGS stress is not as 
crucial as the modelling of the Reynolds stress term in RAS. 

It should be noted here that our treatment of the closure problem, Eqn. (8), is 
what Mason [18] refers to as the Lilly- Deardorff approach. This approach allows us to 
define the filtering operator implicitly through the modelling of the SGS stress term. 
Another approach, referred to as the Leonard approach, is to substitute 



into UaUp to get 

Note that the first term on the right hand side of the above equation can be exphcitly 
calculated if the filter used is itself explicitly defined. Using the above, the SGS stress 
term in Eqn. (8) becomes 



The advantage of the Leonard approach is that it breaks the SGS stress into three 
terms that have physical interpretations in the fiow. The first term on the right 
hand side of Eqn. (9) represents the interaction of large eddies that produce small 
eddy effects and is called the Leonard stress. The second term, called the cross term, 
represents the interaction between the large and small eddies and the third term 



Ma = Ma + 



'a 



(9) 
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represents the interaction of small eddies to produce large eddy effects. It is this 
term that produces the transfer of energy from the small to the large eddies and 
soe is known as the backscatter term. Note that the second term may also produce 
backscatter. The disadvantage of the Leonard apprach is that the filter must be 
explicitly known for the calculation of the Leonard stress term. This may not be a 
trivial task, especially in grid base methods where the filter is often implicitly defined. 
As pointed out by Mason [18] and Ferziger [4], both the Lilly-Deardorff and Leonard 
approaches have been used successfully by researchers. The approach taken for our 
work is that of Lilly-Deardorff. 

2.4 Subgrid Scale Models 

As discussed in the previous section, in filtering the Navier-Stokes equations a new 
term is introduced which leads to a closure problem. To deal with this problem, the 
filtered Navier-Stokes equations are rewritten into the form given in Eqn. (8) and 
the new term. Tap, is then modelled. In this section, we will review two of the most 
commonly used SGS models, namely the Smagorinsky and the Dynamic Subgrid Scale 
models, and introduce our variant of the Smagorinsky model, which we will refer to 
as the Modified Smagorinsky model. 
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2.4.1 Smagorinksy Model 

The Smagorinsky [25] subgrid scale model, 

- ^r,,5„, = -ur + = (10) 

where is called the eddy viscosity (to be derived below) and 

_ 1 / du^ _^ dup\ 
2 \dXj3 dxaj ' 

is known as an eddy viscosity type model because it assumes that the small eddies 
remove energy from the flow through a dissipative process. This precludes the pos- 
sibhty of backscatter (the reintroduction of energy from the small to large eddies) 
which is considered one of the model's weak points. The following derivation of the 
model is that of Ferziger [4] . 

According to turbulence theory energy is introduced at the largest scales and is 
successively transferred to the smaller scales until viscous damping becomes the dom- 
inant effect. In the region where the viscous effect becomes dominant, the turbulent 
energy of the flow is damped out by the transfer of the kinetic energy to internal en- 
ergy. Between these two scales is a region known as the inertial subrange, where there 
is neither signiflcant production nor dissipation of energy. In this region, only the 
inviscid mechanisms are active and so assuming that the transfer of energy is always 
from the large to small scales, the term responsible for the transfer of this energy is 
then the advection term of the Navier-Stokes equations. This allows us to estimate 
the rate of energy transfer to the small scales as the magnitude of the contribution 
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of the advection term to the kinetic energy equation, which is 



1 d 



2dxp 



Since energy is introduced at the largest scales, we estimate the dissipation rate, e, 
as 



where C/l is the velocity scale of the large eddies and L is the integral length scale 
of the turbulent flow. Furthermore, if we assume that the largest subgrid scales are 
much larger than the viscous scales, then 



where Uses is the velocity scale of the small eddies and A is the length of the largest 
subgrid scale eddies, which is also the length scale associated with the filter. 

Now, as stated before, turbulence theory states that energy is transferred from 
the large to small scales. We assume that this process is dissipative in nature with 
respect to the large eddies (i. e. once the energy is lost by the large eddies to the 
small, it cannot be recovered.) So, the SGS model represents this energy transfer as 
effective viscous dissipation. Since it is the smallest resolved scales (of size A) that 
are most influenced by the SGS model, we estimate the effective dissipation as 



e « Ul/L 



(11) 




(12) 



(13) 



Using (12) and (13), we get 



ut oc IMJsGS 
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and using the above together with Eqns. (11) and (12) 

Using 

we get 

ut = ClAlLi\S\ 

where Cs is the Smagorinsky parameter which is introduced to produce the equahty. 
Because of the difficulty in calculating the integral length scale, L, most researchers 
make the simplifying substitution 

which leads to the more common form of the eddy viscosity in the Smagorinsky model, 

UT^{CsAf\S\. (14) 

In our calculations we use as our length scale 

A = (AiA2A3)5 

where Aq, is the grid element length in the appropriate direction and set 

Cs = 0.065. 

Although the Smagorinsky model is the most widely recognized and used SGS 
model, it is not without its problems. One significant problem is that in non- 
homogeneous wall bounded flows, near the surface of the walls, the model does not 
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damp out the eddy viscosity enough to allow kinematic viscosity effects to become 
dominant. In fact, as the wall and hence the viscous sublayer is approached, our 
approximation, (12), in the derivation of the model itself becomes questionable since 
the grid points, hence A, will be in or near the viscous sublayer region. 

Another problem as mentioned before, is that since the Smagorinsky model as- 
sumes that the transfer of energy from large to small scales is a diffusive process, 
Eqn. (13), no mechanism exists to transfer energy from the small to the large scales, 
i. e. there is no backscatter. To show this, we follow the work of Piomelli et al. [22] 
and observe that in the resolved energy transport equation, 

_^ dig'^up) _ d / 2 2- _^ 1 \ 2 dua dua _^ ^ 

dt dxp dxp \ '^aT'ap Jiedx^J Redxi^dxjs ^"'^ 

where — UaUa, one-half of the last term, 

^SGS — TafsSa/S 

which is referred to as the "subgrid-scale dissipation", represents the transfer of 
energy between the large and small scales. Note that 

esGS < 

indicates that the resolved energy decays with respect to the interaction between the 
large and small scales (forward scatter), while 

^SGS > 

means that the resolved energy grows with respect to the interaction of the large 
and small scales (backscatter). In the case of the Smagorinsky model. Equations 
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(10) and (14) cleary show that escs < 0, and so there is no backscatter when the 
Smagorinsky model is used. In their study using DNS data, PiomeUi et al. [22] found 
that for channel flow using a Sharp-Cutoff filter, nearly 50% of the points experienced 
backscatter. For the Gaussian filter, the percentage dropped to about 30% and the 
Top-Hat filter fell inbetween the cutoff and the Gaussian. These results indicate 
that although in the mean, the net effect of the SGS dissipation is forward scatter, 
there is significant backscatter. The importance of backscatter in the SGS model is 
not yet clear. In fully developed channel flows for example, purely dissipative SGS 
stress models have been used successfully. Since forward scatter is dominant in the 
mean, this would seem reasonable. As PiomeUi et al. [22] speculated, it would seem 
reasonable that backscatter effects would be more important for nonequilibrium flows. 

Finally, there is the problem of setting the Smagorinsky parameter, Cs- The 
values used range from about 0.2 for isotropic turbulence to 0.065 for the channel 
flow. Other than for isotropic turbulence, the author is not aware of any systematic 
analysis to determine the parameter. 

2.4.2 Dynamic Subgrid Scale Model 

To overcome the problems of the Smagorinsky model, Germano et al. [6] introduced 
a method of calculating the Smagorinsky parameter interactively during run time 
rather than setting the parameter to a fixed constant at the beginning of the numercial 
calculations. Here, we introduce a variant of Germano's dynamic procedure developed 
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by Piomelli and Liu [23]. The procedure involves applying a second filter, which is 
called the test filter ~, to the filtered momentum equation, Eqn. (7). This results in 
dUa d df d%B 1 92 ~ 

+ ^i^aUp) = ^ h S-^Ua 

at oxp ox a ox 13 Re axp 

where 

Now we define the resolved turbulent stress as 

which represents the contribution of the small resolved scales to the Reynolds stresses. 
Using the Smagorinsky model as our SGS model, we have at the filter and test filter 
level. 

Tap - Ir^-ySaP ~ -2C s \S \S = -2CsB^p (15) 

Tap - \tii5^p ~ -2CsA^\E{S^p = -2CsA^p (16) 

where we have used Cs rather than C| in the eddy viscosity term, Eqn. (14), and A 
and A represents the width of the first filter and the test filter, respectively. Now, 
the SGS stresses at the filter and test filter level are related to the resolved turbulent 
stress through an algebraic identity known as Germano's identity, 

— — "Tap- 

Note that if we use the models for the stresses, Eqns. (15) and (16), then Germano's 
identity is no longer satisfied exactly. Hence, we obtain the residual equation 

EaP = ^af3 + 2CsAai3 — 2CsBaj3- (17) 
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Piomelli's method is to minimize the square of the residual assuming that the Smagorin- 
sky parameter, Cs, is known (denoted by Cg) in the third term on the right side of 
Eqn. (17). So the residual equation becomes 

and minimizing the square of the residual yields 

2 A^ri-^^ri 

The Smagorinsky parameter, Cg is then approximated using a Taylor series expansion 
oiCs. 

One of the advantages of the DSGS model is that there are no parameters to 
determine before the numerical run. Also, the DSGS model is able to reduce the 
eddy viscosity in the laminar regions of the flow. Finally, the DSGS model can allow 
for backscatter in the flow. 

Although Piomelli's version of the DSGS model was formulated to answer some 
of the problems with Germano's original dynamic procedure, there are still some 
unresolved issues with the model. One of the problems is that in the derivation of 
the model, the Smagorinsky parameter was assumed to be independent of the filter 
width (see Eqns. (15) and (16).) However, using DNS data from isotropic turbulence, 
Meneveau [20] found that the instantaneous values of the Smagorinsky parameter 
varied significantly when computed at two different filter widths. 

Another problem with the DSGS model is that, as formulated, it allows Cs and 
hence the eddy viscosity to be negative. Although this allows the model to account 
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for backscatter, if the eddy viscosity remains negative long enough, it can lead to 
numerical instabilities. Researchers have found this to be case and some type of an 
artifical constraint is usually required to prevent this instability. In our case, we use 
the constraint 

where Cmin = —0.01. 

Also, during the execution of the DSGS model, it was found that if we calculate 
Cs at each time step, then Eqn. (18) can and does become unstable. So, to stabilize 
the calculations, we use time averaged values to calculate Cs- The time interval used 
is kept small (100 time steps) to preserve the dynamic nature of the model in time 
as well as in space. 

Finally, although not so much a problem as it is a question, is the issue of the 
double filtering. In the derivation of the DSGS model, it is assumed that filtering 
twice is equivalent to filtering once with some filter. This assumption is used to apply 
the Smagorinsky SGS model to the SGS stress term at the test filter level, Eqn. (16). 
In Appendix B, we derive explicitly the filters that are equivalent to filtering twice 
with the Gaussian and the Top-Hat filter. 

2.4.3 Modified Smagorinsky Model 

One of the problems of the Smagorinsky model as mentioned above, is that the eddy 
viscosity does not damp out as the walls are approached. To overcome this problem 
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we introduce a damping function in the calculation of the Smagorinsky parameter. 
The damping function is constructed under the constraint that the distance from the 
wall cannot be used directly by the model and that the new model cannot significantly 
increase resource usage (e. g. cpu time and memory.) 

To construct our damping function, we first define the instantaneous kinetic energy 

as 

_ 1_ _ 

the instantaneous dissipation rate as 

and the instantaneous Reynolds number as 

ve 

Note that near the wall, if we Taylor expand Ua assuming that the wall location is 
at X = 0, then the only non-zero derivative of the velocity vector will be in the wall 
normal direction, y, assuming the wall to be in the x-z plane. So 



+ ?/■ 



x=0 



+ 

x=0 



dy 

By the no slip boundary condition, TZa(O) = 0, and so 

0{k) = y\ 

Using a similar argument, we have that 

0(e) = 1. 
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So 

Reinst ~ y^- 

Hence we use the instantaneous Reynolds number, Rcinst, to detect the wall. Since 
the Smagorinsky model works well away from the walls, our modified Smagorinsky 
model is 

= {C'sAf\S\ (19) 



where 



C's — C^s < 



for < Reinst < Re" 
1 for Reinst > Re° 



and 



Re l^\_Reins^max- 

Note that 7 and are just model parameters and that [Reins^max represents the 
largest value of Reinst in the compuational domain at a given instant in time. 

As will be shown in the results section, our model is able to detect the walls and 
damp out the eddy viscosity. However, there is nothing in the derivation of the model 
that indicates that Reinst < Re° only near the wall surfaces. In fact, as will be seen 
in the results for the channel and backward facing step problems, Reinst < Re° occurs 
quite often in elements far from the walls. 
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3 Channel Flow 

We use the channel flow as a benchmark to test the three subgrid scale models, because 
of the extensive amount of data, both experimental and theoretical, that is available 
for this type of flow. Since fully developed turbulent channel flow is homogeneous in 
both the streamwise and spanwise directions, we use periodic boundary conditions in 
those directions. 

3.1 Computational Parameters 

The computational domain of the channel consists of 53 x 32 x 16 elements and its 
dimensions are lOh in the streamwise direction (x), 2h in the spanwise direction (z) 
and 2h in the wall normal direction (y), where h (= 0.5) is the channel half height. 
A sketch of the channel is given in Fig. 1. 
For initial data, we use the following: 

u = f(y) + eu* 
< V ^ ev* 
w — ew*, 

where f{y) is a parabolic function zero at the walls and one at the channel center, 
e = 0.1, and u*, v*, and w* e [—1, 1] (uniform distribution.) In our flow experiments, 
we set the viscosity, u — 10~^. So, using mean initial streamwise center line velocity, 
Uq — 1 as our velocity scale, and the channel height, 2h — 1, as our integral length 
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2h 



y=h 



y=-h 



lOh 




Figure 1: Sketch of the channel with relevant dimensions. For our experiment, we 
use h = 0.5 

scale, we estimate the Kolmogorov length scale as 

(i?eo)"2 = 0.001. 



2hUo\ * 



Hence, near the walls, we grade the grid so that the smallest element will have a 
height of O(10~^). A sketch of the graded grid in the xy-plane face of the channel 
is shown in Fig. 2. Note that the grading in the streamwise direction was done to 
reduce the number of grid points needed in the streamwise direction. It is not a result 
of any physics of the flow, as in the wall normal direction. A uniform grid is used in 
the spanwise direction. 



3 CHANNEL FLOW 



33 




Figure 2: Sketch of the xy-plane face of the graded channel grid. 
3.2 Constant Flow Rate 

In assuming periodic boundary conditions in the streamwise direction, the compu- 
tational grid is in essence following the flow down stream. In order to maintain a 
time dependent flow, an external force must be applied to the flow. In the case of 
the channel, a negative streamwise pressure gradient is used to sustain the flow. To 
calculate this pressure gradient, we follow the procedure described by Tran and Mor- 
choisne [27] in which the pressure gradient is found in such a way that the flow rate 
remains constant. Let 

F, = (F,(i),0,0) 
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be the streamwise pressure gradient which may vary in time but is constant in space. 
Using the centerhne velocity, Uq, and channel half height, h, as our characteristic 
velocity and length, the non-dimensionalized momentum equations are 



dUa d dP 1 d'^Ua 



(20) 



where ReQ = ^. 

Since the flow is assumed to be homogenous in the stream and spanwise directions, 
we define our averaging operator to be 



< / >x,z (y) 



1 



LxLz Jo Jo 



f{x, y, z) dxdz, 



(21) 



where L^ and L^ are the streamwise and spanwise lengths of the channel, and the 
mean flow rate as 

dy. (22) 
So, integrating the streamwise momentum equation, Eqn. (20), we get 



1 



h L^Lz J -I Jo Jo 



1 fL. fL^ dU d , dP 



dt dx 



f3 



_ 1 d^U 
dx Rcq dxpdxjj 
c 



E 



D 



> dx. (23) 



We simplify each term in the above equation as follows: 



Using Eqn. (21) and (22), A is simply ^. 



For term B, we break up the sum and simplify each term seperately. The first 



term is 



-iJo Jo ox J -I Jo 



x=0 
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since the flow is periodic in the x direction. Simiharly, since the flow is also 
periodic in the z direction, we have that 

"1 fLz fLx Q 



/I pl^z pl^x r) 
Jo /„ S(^^''^'' = °- 



For the third term, 



/I rLz pLx Q rLx rL^ 

-1 Jo Jo dx^ ^ Jo Jo 



y=l 



= 



y=-l 



by the no-slip boundary condition along the walls. So, B — 0. 

• C must also be zero, for if the mean pressure gradient was not zero, then there 
would be no need to add an external force, F^, to drive the flow. 

• For D, we break up the sum as we did for B. So 



I If 

J-iJo Jo 







x=0 



dx^ J-i Jo dx 

since the flow is periodic in the x direction. A similiar argument can be used 
to show that 



1 rLz rL^ g^U 



0. 



-1 ^0 ^0 dz'^ 
For the third term, we have that 

L, 1 /•! f^^ f^- 1 d^U L^ 1 r^d'^<U>^, 



h L^L^ J— 



1 fi^z n 
1 Jo Jo 



h Re, 



-J 

-^-1 



dy 



L, 1 d<U>, 



h Rcq dy 



y=l 



y=-l 



Since the mean channel flow is symmetric about the channel center line, y — 0, 



d<U>, 
dy 



y=l 



d<U>, 
dy 



y=-l 
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So, 



^ L, 1 d<U >:c,z 



h Reo dy 



h Reo dy 
Finally, since F^ is constant in space, E simplifies to 



y=-l 



h L^Lz J-i Jo Jo ^ h 



Using the above, Eqn. (23) can be rewritten as 



dQ ^L^ 2 



where ul — 



1 d<U> 



X.Z 



is the mean wall shear stress. So, as the flow changes 

y=-l 



"t Reo dy 

from a laminar state to a turbulent state, increases and if the pressure gradient 
is kept constant at the laminar value, Eqn. (24) says that Q will decrease. To main- 
tain a constant flow rate, it is then necessary to balance the pressure gradient to 
the mean wall shear stress. However, to calculate the mean wall shear stress would 
require a large grid to obtain the mean values in the stream and spanwise directions. 
This makes balancing the pressure gradient with the mean wall shear stress compu- 
tationally too expensive. So, F^ is found by relating the pressure gradient to the flow 
rate fluctuations from the constant initial value Qo, 

where ai — 2At, a2 — At, and n and n + 1 represent the old and new time states. 
Note that the formula given in the reference, [27], contains typographical errors. 
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3.3 Numerical Results 

Using the mean initial centerline velocity as our characteristic velocity and the channel 
half height as our characteristic length scale, the Reynolds number for our numerical 
experiments is Re^ — 5000. Note that in the following discussion that all mean 
calculations are based on averaging in the homogeneous directions (x and z) as well as 
time. Also, note that we use as our characteristic velocity for non-dimensionalization, 
the wall shear velocity 



where is the wall shear stress and use as our non-dimensional wall distance 



Fig. 3 is a plot of the total kinetic energy using the three different models. As 
can be seen in the figure, by simulation time t — 150, the flow has stabilized enough 
to begin the time averaging calculations. Note that in the DSGS model run, the 
Smagorinsky model is used while the flow is developing until time t — 75. This is the 
reason for the initial overlap of energy between the Smagorinsky and DSGS results 
in Fig. 3. 

One of the problems with the Smagorinsky model discussed above is that the eddy 
viscosity is not damped as the wall and the viscous sublayer is approached. This can 
clearly be seen in Fig. 4, the sketch of the time averaged eddy viscosity along the hues 
X = 5h and z = 1.0625/i, and x = 9h and z = 1.0625/i. Note that along these lines. 
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Figure 3: Plot of the total kinetic engery for the channel flow using the Smagorinsky, 
DSGS, and modified Smagorinsky subgrid scale models. At time t = 150, a running 
total of the fluid variables are started to calculate the various mean properties of the 
flow. 

the magnatude of the eddy viscosity is approximately 75% of the kinematic viscosity 
at the walls, when using the Smagorinsky model. The plots also show that both 
the DSGS and the modified Smagorinsky model damp out the eddy viscosity as the 
wall is approached. Fig. 4 also indicates that the damping function in the modified 
Smagorinsky model not only damps out the eddy viscosity near the walls but also 
causes some damping far from the walls. However, the plots show that the damping 
along the central region of the channel is not as great as the damping near the walls. 
Fig. 5, a snap shot of the eddy viscosity along the line x — 5h and z — 1.0625/i at 
simulation times t — 150 and t — 300, shows that at any given time, the eddy viscosity 
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can be over a magnatude greater than the kinematic viscosity when using the DSGS 
model. Furthermore, Fig. 5b shows that at any given time, the eddy viscosity can be 
negative (which can be interpreted as backscatter) when using the DSGS model. 

Fig. 6, a log-linear plot of the non-dimensionalized mean velocity, indicates that 
the modified Smagorinsky data is much closer to the logarithmic friction law than the 
Smagorinsky model. Note that we use as our additive constant in the log law, 5.5, as 
in Kim, Moin, and Moser's work [12] since our Reynolds number, based on the wall 
shear velocity, is Rct ~ 180 for our channel computations. 

Figures 7-9 are plots of the root-mean-square of the velocity fluctuations. Note 
that the DNS results to which we compare our data is that of Kim et al. [12] and the 
experimental data is that of Kreplin and Eckelmann [13] . It is not very surprising that 
the velocity fluctuation properties of our LES computations do not match the DNS 
and experimental results exactly, because of the way in which the fluctuations are 
computed. If we let < • >t,x,z be the average in time and the homogeneous directions 
and ' denotes the statistical fluctuations with zero mean, then note that 

{u'af =< Ua >t%,^ -<ul> . 

However, for our LES cases, the only values we have are at the filter level, Ua- Hence 
the only fluctuation we can calculate is 

(Kf =< "^a >t,:c,z -<^l >t,x,z 

where u'^ represents the statistical fluctuations in the flltered fleld. Note that u'^ is 
well defined (in the sense that it will not be identically zero) since, as discussed in the 



3 CHANNEL FLOW 40 

filter section, only the Sharp-Cutoff filter has compact support in wave space and so 
for filters such as the Top-Hat filter, Ua will still contain some small scale motions. 
Also, the grading of the grid near the walls imphes that the filter width decreases 
near the walls and so will contain more of the small scale motions. Hence, we would 
expect u'^ to be closer to u'^ near the walls, were the filter width is much smaller than 
toward the channel center. Figures 7-9 show that this is indeed the case. 

Fig. 10 is a plot of the Reynolds shear stress normalized by the wall shear velocity. 
As discussed above, we note that the LES results are much closer to the DNS results 
near the wall. Also note that the DSGS and modified Smagorinsky are closer to the 
DNS results than using standard Smagorinsky model. 

Figures 11-13 are plots of the root- mean-square vorticity fiuctuations normalized 
by the mean wall shear. Since vorticity is mostly associated with small scale motions, 
it is not surprising that the LES vorticity data are all lower than the DNS data for 
the reason given above. Note that near the wall, at least one of the LES data is 
close to the DNS results. Unfortunately, there is no consistency as to which model 
is closest to the DNS results. For x vorticity fiuctuations, the Smagorinsky model is 
closest to the DNS data near the wall, the modified Smagorinsky is the furthest off 
from the DNS data, and the DSGS model is between the other two. However, for the 
z vorticity, note that the order of performance of the models changes to the DSGS 
beging the closest, followed by the modified Smagorinsky and then the Smagorinsky. 
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Figure 4: Time averaged plots of the eddy viscosity normalized by the kinematic vis- 
cosity, (a) is a plot along the line x — 5h and z — 1.0625/i. (b) is a plot along the 
line X — 9h and z — 1.0625/i. 




Figure 5: Instantaneous snap shot of the eddy viscosity normalized by the kinematic 
viscosity along the line x = 5h and z = 1.0625/i. (a) is using the Smagorinsky model, 
(h) is the DSGS model, and (c) is the modified Smagorinsky model. 
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Figure 6: Log-linear plot of the mean velocity normalized by the wall shear velocity, 
Ur, for the channel flow using the Smagorinsky, DSGS, and modified Smagorinsky 
subgrid scale models. 




Figure 7: Root-mean-square of the streamwise velocity fluctuations normalized by the 
wall shear velocity. Note that dns data is from Kim et al. [12] and the experimental 
data, exp, is from Kreplin and Eckelmann [13]. 
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Figure 8: Root-mean- square of the wall normal velocity fluctuations normalized by the 
wall shear velocity. Note that dns data is from Kim et al. [12] and the experimental 
data, exp, is from Kreplin and Eckelmann [13]. 




Figure 9: Root-mean-square of the spanmwise velocity fluctuations normalized by the 
wall shear velocity. Note that dns data is from Kim et al. [12] and the experimental 
data, exp, is from Kreplin and Eckelmann [13]. 
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Figure 10: Reynolds shear stress normalized by the wall shear velocity. Note that dns 
data is from Kim et al. [12] and the experimental data, exp, is from Eckelmann [3] 
at two different Reynolds numbers, Rcr — 142 and Rcr — 208. 
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Figure 11: Root-mean- square vorticity (x) fluctuations normalized by the mean wall 
shear. Note that dns data is from Kim et al. [12]. 
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Figure 12: Root-mean-square vorticity (y) fluctuations normalized by the mean wall 
shear. Note that dns data is from Kim et al. [12]. 
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Figure 13: Root-mean-square vorticity (z) fluctuations normalized by the mean wall 
shear. Note that dns data is from Kim et al. [12]. 
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4 Flow over a Backward Facing Step 

To study the effectiveness of the various SGS models on turbulent flows with sep- 
aration and reattachment, the flow over a backward facing step is used because of 
its geometric simplicity. In our modeling of this flow, we use periodic boundary 
conditions in the spanwise direction. 

4.1 Computational Parameters 

The computational domain of our backward facing step problem consists of one main 
channel (post expansion section) with the step located at one end of the channel, see 
Fig. 14. The discretization of the channel consists of 88 x 32 x 16 elements and the 
dimensions of the channel are approximately 20h in the streamwise direction and 2h 
in the spanwise and wall normal directions, where h = 0.45833 is the step height. 

At the inlet, the fluid is forced by setting the streamwise velocity to be 1 along 
the inlet plane. Thus the Reynolds number, based on the inlet velocity and the step 
height, is Rch = 4583. 

4.2 Numerical Results 

As in the channel flow simulations, the time averaging calculations are started once 
the total kinetic energy has reached a quasi-equihbrium state, see Fig. 15. Figure 16 
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20h 




Figure 14: Sketch of the backward facing step and channel with approximate relevant 
dimensions. For our experiment, the step height is h — 0.45833 

is a plot of the mean skin friction coefficient, 



along the bottom wall, where Uo is the inlet velocity. According to Le & Moin [15], 
the high |C/| in the recirculation region is the result of not having an entry section. 

Figures 17 and 18 are comparisons of the mean streamwise velocity at various 
locations in the recirculation, reattachment, and recovery regions. In the reattach- 
ment and recovery regions, the data indicate that near the wall all the LES models 
overpredicted the mean streamwise velocity. In the region of the secondary bubble, 
the skin friction data (Fig. 16) and the mean velocity data (Fig. 17) indicate that for 
the Smagorinsky model and the modified Smagorinsky model, the secondary bubble 
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is not as well defined as in the DNS and DSGS cases. 

Figures 19-22 are the profiles of the various Reynolds stress components. Figure 
19 indicates that in the reattachment region, all the SGS models predicted a greater 
streamwise turbulence intensity near the wall than the DNS results. According to 
Le and Moin [15], their results were lower than the experimental data of Jovic and 
Driver [11]. Also, the data from AkselvoU and Moin [1], indicate that their LES also 
predicted higher turbulence intensities, although not as high as our results. In the 
recirculation region, Le and Moin [15] reported that the maximum of all Reynolds 
stress components occurred around the step height. As can be seen in Figures 19-22, 
the maximum of all Reynolds stress components in the LES cases occur before the 
step height and does not peak as sharply as in the DNS case. 
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Figure 15: Plot of the total kinetic energy for the backward facing step using the 
Smagorinsky, DSGS, and modified Smagorinsky subgrid scale models. At time t — 
150; a running total of the fluid variables is started to calculate the various mean 
properties of the flow. 
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Figure 16: Skin friction coefficient along the bottom wall of the post expansion section 
of the backward facing step problem, dns: direct numercial simulation data of Le and 
Moin [15], smag.: original Smagorinsky subgrid scale model, dsgs: dynamic subgrid 
scale model, and mod. smag.: modified Smagorinsky model. 
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Figure 17: Mean streamwise velocity profiles in the post expansion section of the 
backward facing step problem, dns: direct numercial simulation data of Le and Main 
[15], smag.: original Smagorinsky subgrid scale model, dsgs: dynamic subgrid scale 
model, and mod. smag.: modified Smagorinsky model. 
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Figure 18: Mean streamwise velocity profiles in the post expansion section of the 
backward facing step problem, dns: direct numercial simulation data of Le and Main 
[15], smag.: original Smagorinsky subgrid scale model, dsgs: dynamic subgrid scale 
model, and mod. smag.: modified Smagorinsky model. 
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Figure 19: Streamwise turbulence intensity profiles in the post expansion section of 
the backward facing step problem. DNS: direct numercial simulation data of he and 
Moin [15], Smagorinsky: original Smagorinsky subgrid scale model, DSGS: dynamic 
subgrid scale model, and Mod. Smagorinsky: modified Smagorinsky model. 
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Figure 20: Wall normal turbulence intensity profiles in the post expansion section of 
the backward facing step problem. DNS: direct numercial simulation data of he and 
Moin [15], Smagorinsky: original Smagorinsky subgrid scale model, DSGS: dynamic 
subgrid scale model, and Mod. Smagorinsky: modified Smagorinsky model. 
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Figure 21: Spanwise turbulence intensity profiles in the post expansion section of the 
backward facing step problem. DNS: direct numercial simulation data of Le and Moin 
[15], Smagorinsky: original Smagorinsky subgrid scale model, DSGS: dynamic subgrid 
scale model, and Mod. Smagorinsky: modified Smagorinsky model. 
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Figure 22: Reynolds shear stress profiles in the post expansion section of the backward 
facing step problem. DNS: direct numercial simulation data of Le and Moin [15], 
Smagorinsky: original Smagorinsky subgrid scale model, DSGS: dynamic subgrid scale 
model, and Mod. Smagorinsky: modified Smagorinsky model. 
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5 Computer Software 

The numerical package used in this study is known as Hydra. This software was 
developed by Mark Christon for Lawrence Livermore National Laboratory. The finite 
element code was in turn developed on the works of Gresho et al. [9], [10], [8]. The 
original Smagorinsky turbulence model was incorporated into the code based on the 
work of Rose McCallen. Also a disscussion of the motivation for using the finite 
element method is also presented by McCallen [19]. 

Although Hydra has several options for solving the Pressure Poisson Equation, 
it was found that for our turbulent flow cases the only solver that worked within a 
reasonable period of time was a version of the direct solver. This together with our 
limited computer resources restricted our mesh size to the 88 x 32 x 16 grid that was 
used for the above results. 
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6 Conclusions 

One of the unique characteristics of LES as compared to other methods of deahng 
with turbulent flows is the idea of filtering. As was mentioned before, this concept of 
filtering introduces a new error to the system which is the result of the noncommu- 
tative nature of the filtering operation and differentiation in graded computational 
grids. Most of the focus in studying this error has been on the filter kernal function 
G and trying to obtain a G that would allow control over the order of the interchange 
error. As we used FEM as our numerical method, we had knowledge of the general 
shape of the solution (element wise) and so decided to focus on trying to use this ex- 
tra information. We were able to show that under certain conditions the interchange 
error is second order. 

As we have noted previously, one of the weaknesses of the Smagorinsky subgrid 
scale model is its inability to damp out the eddy viscosity near walls. The purpose 
of this project was to develope and implement a new subgrid scale model that would 
address this problem. To test the new model two turbulent flows where studied, the 
channel flow and the backward facing step. The channel flow was chosen because of 
the extensive amount of data, experimental, computational, and analytical, that is 
available for that type of flow. The backward facing step was chosen because of its geo- 
metric simplicity for a flow with separation and reattachment. Also, the Smagorinsky 
and a version of the Dynamic Subgird Scale model were also implemented for com- 
parison. These two models were chosen for comparison for two reasons. First, they 
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represent the two SGS models that are in common use in LES and secondly, their 
treatment of the eddy viscosity near walls are very different. As mentioned above, the 
Smagorinsky model does not damp out the eddy viscosity near the walls whereas the 
DSGS model does damp them out. Although the DSGS model does damp out eddy 
viscosity near the walls, some problems with the model, as mentioned previously, are 
that it requires double filtering and the eddy viscosity becomes negative and leads to 
numerical instability if left unchecked. 

The results of our computations with the modified Smagorinsky model were com- 
parable to the DSGS model computations and generally performed better than the 
original Smagorinsky model. Our eddy viscosity data from the channel flow indicate 
that our model and the DSGS were indeed able to detect the walls of the channel and 
reduce the eddy viscosity accordingly. 

The next step in the continuation of this work should be to develope a better 
numerical solver that would allow a finer grid and, in the case of the backward facing 
step, allow an inlet channel to also be used. 
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A Filter /Derivative Interchange Error 

In this section we analyze the error in the interchange of filtering and differentiation 
by numerically calculating the filter of the derivative and the derivative of the filter 
and comparing the difference between the two. We begin by first creating our graded 
mesh. The mesh is generated by taking an interval [a,b], in our case [0,1], and 
partioning it into n even intervals. The graded mesh is then generated by using the 
following transformation: 



where h — - and a is some constant, which in our case is a = 5. If we define the 
filter width, A (a;), to be twice the length of the smallest element containing x, then 
we have 



.a{a+ih) 



- 1 



A{xi) 




- 1 - 1 

,a{a+ih)f-\ p—ah\' 



6°" - I 



- 1 




So 



1 - e 



,—ah 



A{x) = 2 x{l- e-'^^) + 



- 1 



(25) 



and 



A'(a;) = 2(1-6-"'*). 



(26) 
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Taking the Taylor expansion of the exponential function about in Eqn. (25) and 
(26), it is clear that the order of the filter width and its derivative is h, i. e. 0{A{x)) — 
0{A'{x)) = 0{h). 

The numerical differentiation is performed using a fourth order method and the 
integration is carried out using the Composite Simpson's Rule over each element in 
the filter support. This ensures that the numerical integration is also at least fourth 
order in h. We use a truncated Gaussian as our filter and two different filter widths 
were tested, Ai(a;i) = 2{xi — Xi^i) and A2(a;j) = 2{xi — Xi_2)- 

To test the validity of our assertion that the numerical differentiation and inte- 
gration methods used are in fact fourth order, we observe the following. In section 
2.2, page 11 we found the interchange error to be 

A\x) [\yG{yH{x-yA{x))dy. . 

Note that if our original function is a line, e. g. u{x) — lOx + 1, then our polynomial 
interpolations will all be exact, i. e. Ue{x) — u{x) everywhere. Furthermore, the 
derivative is constant and the same in every element, e. g. for the above linear u, 
u'^{x) — 10. So the interchange error should be 0. Hence the only error observed 
with the linear u will be the error in the calculation of the derivatives and integrals. 
Table 1 indicates fourth order behavior and so our claim of fourth order accuracy of 
the numerical differentiation and integration is justified. 

We test the interchange error by using a sinusoidal source function 

u{x) = sin(27ra;) + 1. 



du . , du . , 

Tx^""^ - Tx^""^ 
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Table 2 clearly indicates that the interchage error is second order. Hence our claim 
that the interchange error is second order if 0{A'{x)) — 0{A{x)) is justified. 
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f (^) - ^ 




linear interpolation 


quadratic interpolation 


h 






Ai(a;) 


A2{x) 


0.05 


1.20e-03 


1.20e-03 


1.20e-03 


1.20e-03 


0.005 


1.19e-07 


1.06e-07 


1.19e-07 


1.06e-07 


0.0005 


1.72e-10 


1.68e-10 


1.70e-10 


1.68e-10 



Table 1: Max norm of the diff. /filter interchange error, linear source function u{x) = 
lOx + 1 









linear interpolation 


quadratic interpolation 
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Ai(a;) 


A^ix) 


Ai(a;) 


A^ix) 


0.01 


1.07e-02 


4.31e-02 


1.37e-02 


4.94e-02 


0.001 


1.20e-04 


4.97e-04 


1.55e-04 


6.13e-04 


0.0001 


1.60e-06 


5.00e-06 


1.56e-06 


6.18e-06 



Table 2: Max norm of the diff. /filter interchange error, sine source function u{x) — 
sin(27ra;) + 1 
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B The Effect of Double Filtering 



The following is a verification of the assumption, in Germano's derivation of his 
Dynamic Subgrid Scale Model, that filtering twice is equivalent to filtering once with 
some filter. In this paper we look at two filters that are often mentioned with regard 
to LES. Namely the Gaussian filter and the Top-Hat filter. 

B.l Gaussian Filter 

Given a Gaussian function 



which is a more general case of the familiar Gaussian filter in LES. It will be easier 
to insert the parameters later to change this into the proper Gaussian filter. Let 




the Fourier Transform of G{t) denoted as G{k) is 




Now we look at the following integral. 






and 



G2{t) 
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then we define 



u[x) 



,) = j - i)di = u{x) * 

and 

u{x) — u{x) * G2{x) 

So the Time Convolution Theorem, which states that the Fourier Transform of the 
convolution of two functions is equal to the product of the Fourier Transform of each 
of the functions, gives us 

u{k) ^ u{k)Gx{k) 

and 

Ji{k) = u{k)G2{k) 

= u{k)Gi{k)G2{k) 



Hence using the inverse transform, we get 



u{k)e 



~ TT 1 /■ -fc^ ( ^ \ ^\ 

u{x) = -^=— u{k)e ' y'-' '^^Je^'^'^dk 
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Now if we let 
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so 



and 



then 



So 



as 



aia2 



Gs{k) = \ —e 
V c^s 



= ^ ^ f u{0 [ G,{k)e^'^^-^Ukd^ 



Now the Gaussian filter for LES is 



So if we let 



6-2 (x) = 



'6 1 -6^ 
TT A^^ 



TT A2 
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then CKi — a2 — and 0:3 = 77^7-^ • Hence if we go back and let • and • be the 



proper Gaussian filters defined with the above kernals, we have that 



u{x) — [u{x) * Gi{x)] * G2{x) 
6 0r A1A2 



TT A1A2 ^6(Af + A 
/ ^(0^ 



e 



^ ^Af + Ai 



where Gs{x) is the Gaussian filter with filter width A = WAf + A^. 



B.2 Top-Hat Filter 

The Top-Hat filter is defined as 

Pa (x) 



1_ 

A 



1 + — 



A 



1 - 



A 





\x\ < 



: otherwise 
We also define a triangle filter as follows. 



-A<x <0 
0<x < A 
otherwise 



Now the Fourier Transform of the above functions are 

Pa(/c) = 



2sin(f) 



and 



kA 



4sin^(^) 
A^k^ 
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If we define • to be tfie Top-Hat filter and • to be the triangle filter, then note that if 
we filter twice with the same (width) Top-Hat filter we get 

u(x) — u* va * Pa 

2 2 

So the Fourier Transform of ^ is 



u{k) — u{k)p A{k)p A(k) 



u{k)- 



4 sin 



2/fcA^ 



Hence 



u{x) — 



I !'mk)e"'^dk 



Jm 
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2 / fc A ^ 



2 J ^ik(x-^) 
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4sin2(f ) 
— — ( 



1 r4sin^(^) 



u(x] 



where 



1 + ^ : X - A < C < X 



A 



1-^ : X < ^ < x + A 
: otherwise 

So filtering with the same Top-Hat filter twice is equivalent to filtering once with 
the triangle filter, q/^. 
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